Critical slowing down in polynomial time algorithms 



A. Alan Middleton 
Department of Physics, Syracuse University, Syracuse, NY 13244 
(February 1, 2008) 

Combinatorial optimization algorithms which compute exact ground state configurations in dis- 
ordered magnets are seen to exhibit critical slowing down at zero temperature phase transitions. 
Using arguments based on the physical picture of the model, including vanishing stiffness on scales 
beyond the correlation length and the ground state degeneracy, the number of operations carried 
out by one such algorithm, the push-relabel algorithm for the random field Ising model, can be 
estimated. Some scaling can also be predicted for the 2D spin glass. 



There has long been a close link between the concepts 
of statistical physics and the algorithms used to simu- 
late condensed matter systems. The fundamental con- 
nection is the mathematics of graphs, which is applied 
analytically, for example, to compute high temperature 
series. Using sophisticated connectivity algorithms bor- 
rowed from computer science, one can numerically calcu- 
late quantities in percolation to a high precision |^ . The 
Fortuin-Kastelyn cluster representation for magnets 
is the basis for the Swendsen-Wang algorithm By 
implementing nonlocal dynamics, these algorithms can 
greatly reduce the time needed for simulations near a 
phase transition. The study of disordered systems, such 
as spin glasses, pinned vortices in superconductors, and 
random field magnets, lead to the introduction of graphs 
with nonuniform edges Early in the study of dis- 

ordered systems, it was realized that the study of such 
graphs is directly related to issues of computational com- 
plexity. In some cases |^ , computing the ground state of 
a disordered material is computationally intractable, as 
the relevant optimization questions on a graph are NP- 
hard Quite interestingly, some computationally in- 
tractable problems have been found to have phase tran- 
sitions j^. For example, given an ensemble of logical 
expressions characterized by the number of Boolean vari- 
ables N and clauses M, the fraction that are satisfiable 
can exhibit finite size scaling about a critical value of 
M/N. 

In this letter, I present results on the behavior of 
ground state algorithms near phase transitions in two 
models, the random field Ising magnet (RFIM) and the 
2D spin glass (2DSG), which, in contrast with problems 
such as satisfiability, are always solvable in time polyno- 
mial in the size of the graphs. These phase transitions 
lead to singularities in the mean time to solve the graph 
optimization problems. For the RFIM, a close connection 
is made between the critical slowing down of the ground 
state algorithm, the correlation lengths, and the degen- 
eracy of the ground state in the thermodynamic limit. 
Numerical evidence is presented that suggests that the 
dynamic critical exponent is z ~ 1 for the RFIM; scal- 



ing arguments suggest z > 1. The behavior far from the 
transition can also be explained. The algorithm for the 
2DSG is more difficult to analyze, but the observed uni- 
form time per spin in the ferromagnetic (FM) phase is 
seen to be natural. 

Generally, the dynamics at low temperature T is ex- 
ceedingly slow in disordered magnets, as seen in experi- 
ment and in Monte Carlo simulations using local moves 
[Q. The glassy dynamical behavior is attributed to the 
complex structure of the energy functional. The free en- 
ergy barriers to equilibration at a length scale £ scale 
as , so that the time to equilibrate a portion of the 
sample of size t''' are expected to scale as ~ exp(£'^/r). 
However, for some random magnets there are algorithms 
which take time polynomial in TV to find an exact ground 
state. The process of finding the ground state uses "non- 
physical" configurations or moves: a local search or sim- 
ulated annealing that uses only physical configurations 
and local moves is hindered by the large energy barriers. 

When using local Monte Carlo moves at finite T to 
model uniform magnets, the relaxation or correlation 
time at continuous transitions scales as L^, with z\oc > 
^/v 1^. Nonlocal cluster moves, such as used in the 
Swendsen-Wang algorithm, can reduce the critical expo- 
nent z, with Zc\ > a/v with a and v the exponents 
for the heat capacity and correlation length, respectively. 
As the algorithms used to find the ground states of dis- 
ordered magnets do not use local moves and are not de- 
signed to utilize clusters related to a critical point, it is 
less clear how many operations are required. Polynomial 
bounds on the worst case behavior of these algorithms 
do exist. For graphs with n vertices and m edges, the 
highest level version of the push-relabel (PR) algorithm, 
used here for the RFIM, will use 0{n^y/rn) time. Useful 
algorithms for general matching jl^ , applied here to the 
2DSG, use from 0{n'^^^) to 0{nP) operations, up to log- 
arithmic corrections, assuming to cx n. In practice, how- 
ever, for typical disorder realizations of physical interest 
in finite dimension d, these algorithms are much faster, 
with the average running time for many problems scal- 
ing roughly as iV^, with q typically in the range 1.1 — 1.3 
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1^,0. For systems with a single T = fixed point, such 
as the elastic medium in a random potential [0, this 
scaling for the time would not be expected to vary with 
parameters. 

The T = RFIM has a phase transition that has 
been extensively studied using a mapping of 
the ground state to the optimization problem max-flow 
0. The Gaussian RFIM has Hamiltonian Hr = 
— JJ2<ij> ^i^j ~ Si hiSi, where the allowed spin values 
on N = L'^ lattice sites i are Si = ±1, the ferromagnetic 
couplings J are positive, and the random fields hi are 
Gaussian distributed, with mean and variance A^J^. 
In the paramagnetic phase (A > Ac w 2.27), the spins 
are pinned by the external fields hi and the magnitude of 
the net magnetization m = ^ s; vanishes as L — s- 0. 
In the ferromagnetic phase (A < Ac), the ferromagnetic 
coupling J dominates and |to| has a finite limit, with the 
sign of m determined roughly by the total of ^ hi (some 
spins will be reversed by fluctuations in hi ) Results for 
the time for this algorithm to find the ground state in a 
3D cubic lattice with periodic boundaries are plotted as a 
function of A and linear system size L in Fig. |l]. Plots of 
the number of primitive operations executed are nearly 
identical in form. Near the phase transition separating 
the ferromagnetic from paramagnetic phase, the time to 
find the ground state shows a characteristic critical slow- 
ing down. This notable effect reflects the deep connection 
between the dynamics of the ground-state algorithm and 
the physics of the RFIM. 

The application of max-flow algorithms to the RFIM 
is well established, but to make the connection between 
scaling in the RFIM and algorithm timing, it is useful to 
review the algorithm. The network flow algorithm gen- 
erally used to solve the RFIM, because of its speed, is 
the PR algorithm of Tarjan and Goldberg |jl^. The al- 
gorithm described here uses a modification that removes 
the need for source and sink nodes ||l^ , reducing memory 
usage and also clarifying the physical connection. 

The modified PR algorithm starts by assigning an "ex- 
cess" Xi to each lattice site i, with Xi = hi. Residual 
capacity variables Vij between neighboring sites are ini- 
tially set to J. A height variable Ui is then assigned to 
each node via a global update step. In this global update, 
the value of Ui at each site in the set T — {i\xj < 0} of 
negative excess sites is set to zero. Sites with a;^ > have 
Ui set to the length of the shortest path, via edges with 
positive capacity, from i to T. 

The ground state is found by successively rearranging 
the excesses Xi, via "push" operations, and updating the 
heights, via "relabel" operations. When no more pushes 
or relabels are possible, a final global update determines 
the ground state: those sites which are path connected 
by bonds with Ty > to T have Si = —1, while the 
sites which are disconnected from T have Si = 1. A push 
operation moves excess from a site i to a lower height 



neighbor j, if possible, that is, whenever Xi > 0, > 
and Uj — Ui~l. In a, push, the working variables are mod- 
ified according to Xi Xi — S, Xj — s- Xj + S, rij — > Vij — S, 
and Tji Tji + 6, with S = min{xi,rij). Push opera- 
tions tend to move the positive excess towards sites in 
T. When Xi > and no push is possible, the site is 
relabeled, with Ui increased to 1 -I- max{j|^.^.>o| . In 
addition, if a set of highest sites U become isolated, with 
Ui > Uj + 1, for all i G U and all j ^ U, the height Ui 
for all i Cz U is increased to its maximum value, N, as 
these sites will always be isolated from the negative ex- 
cess nodes. A proof of the correctness of the PR flow 
algorithm can be found in standard textbooks and 
its application to the RFIM is well known Periodic 
global updates, here applied every TV relabels, are often 
crucial to the practical speed of the algorithm. The high- 
est site heuristic is used here, which applies pushes and 
relabels where Ui is maximal and Xi > 0. 

The PR algorithm is intuitively appealing: when the 
initial capacities within a region are large, the excess can 
be rearranged so that the positive excesses cancel the 
negative excesses as much as possible. The remaining 
excess values having the sign of the original total excess 
for the region. As r,y = J and Xi = hi initially, large 
capacities correspond to the ferromagnetic bonds being 
strong enough to favor alignments of the spins, with the 
spin direction given by the sign of the total hi for the re- 
gion. If the initial capacities are not large enough (weaker 
J), the regions align independently, according to the lo- 
cal field, and the excesses do not cancel. The number of 
steps needed to move excess across the diameter of a re- 
gion via push operations is bounded below by the linear 
size of the region. Note that the running time in the FM 
phase is somewhat dependent on the majority magneti- 
zation, due to the up/down asymmetry of the algorithm. 

I now argue that these correspondences can be used to 
bound the running time of the algorithm, using the phys- 
ical properties of the RFIM, particularly the behavior of 
the ground state degeneracy in the thermodynamic limit. 
Recent work jl^ has given strong numerical evidence of 
insensitivity of the interior spins in the ground state to 
boundary conditions, when A > Ac. (This is in contrast 
with the scenario for, e.g., mean field Q| or highly disor- 
dered [TtI spin glasses in c? > 8, where the entire solution 
is sensitive to the volume.) Most importantly, this im- 
plies that the ground state solution is determined by the 
hi within a volume typically of the size of the correlation 
volume However, in the FM phase, the interior is sen- 
sitive to the boundaries and a finite fraction of the spins 
flip infinitely often as the sample volume is increased, as 
the net magnetization is given by a coarse-grained global 
sum of the hi. 

First, consider the case A ^ 2d, where, in the ground 
state, almost all spins satisfy Si = sgn{hi). In the al- 
gorithm, the initial positive excesses at sites i where 
hi > 2d J remain positive, as the total capacity of bonds 
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leaving the site is only 2dJ, so that push operations can 
only reduce hi by 2d J. The sites with excess Xi < —2d J 
also always have negative excess. There will be some 
small rearrangement, but only locally, and the sign of 
the excess will change only at a very few sites during 
the execution of the algorithm, which terminates after a 
number of operations cx N (up to logarithmic corrections 

A similar scenario holds, but at scales ^, for A > Ac. 
The algorithm establishes the boundaries of correlation 
volumes by rearrangement of excess over distances of 
scale ^. Further rearrangement is blocked by the effec- 
tive decrease in residual capacity with scale (as the stiff- 
ness, corresponding to the scale dependent J, decreases 
rapidly on scales greater than ^ The question to 

be answered, then, is how long the algorithm takes, per 
spin, to redistribute excess on scale ^. The number of 
push and relabel operations in a volume is bounded 
below by For the excess to be pushed over a dis- 

tance f , the relative heights must differ by at least ^, so 
that at least 0{(,'^'^'^) relabels R must be performed for 
each correlation volume (global updates lead to height 
changes, but these do not appear to affect the scaling, 
empirically.) This gives the estimate R ^ L^^. This 
scaling is consistent with numerical results, up to loga- 
rithmic corrections, for d = 1,3. The inset to Fig. |l| show 
a scaling plot for R, for example, with {RL^^Y^'^S plot- 
ted as a function of the finite size scaling variable SL^^'^ , 
S = (A-Ac)/Ac, with the values v = 1.37 and A^ = 2.27 
fixed parameters, determined independently . A fit in 
d = 1, with 1^ — 2 and Ac — also fixed by known 
results, describes the data well for L < 5 x 10^, with 
R ^ ln(L/^) (without the global relabeling heuristic.) 

One limit where the asymptotic time for the algo- 
rithm can be described more precisely is when A <^ 
[Lln{L)]^^/^ .) Here, the capacities do not limit the re- 
arrangement of excess: the final state is either Xi > 
everywhere or Xi < everywhere, corresponding to a 
uniform Si — ±1 state, according to whether ^ ft-i > 
or ^ < 0, respectively. The dynamics of the PR algo- 
rithm is set by the fluctuations in the total hi in regions 
at each length scale. At any time scale of the compu- 
tation, the positive excess will be pushed towards the 
nearest negative excess region, with the distribution of 
excess negative or positive over a particular volume with 
sites, as the rearrangement of excess will have been 
completed over shorter lengths at an earlier time scale. 
Generally, but especially given sufficiently frequent global 
updates, the sites with greatest Ui will be furthest from 
the set T. As the highest sites are examined for push- 
ing, the excess will be moved from these sites to the next 
lower level. These sites will then have their excess moved 
to the next lower level and the algorithm will "sweep" the 
excess through the volume of diameter « 2£. This will 
establish a distribution of excess that will have uniform 
sign over a region of size 2£. In this fashion, the algorithm 



will solve for the sum of hi recursively. The number of 
steps at each stage will be L'^ and O(lnL) stages will be 
required, giving a total running time ~ L'^ln(L). This 
result is consistent with the numerical timings for very 
small A, with the data for RL^'^ linear in \n{L) over sam- 
ple dimensions 16 < L < 128 to within the 1% numerical 
error for d = 3 and over 10^ < L < 5 x 10^ to within the 
same error for d = 1. Coarsening of the height variables 
during the algorithm is displayed in Fig. ^ 

To indicate the generality of critical slowing down and 
partial arguments for other systems, it is useful to com- 
pare the RFIM results with results for the 2DSG. The 
Hamiltonian is Hs = J2<ij> JijSiSj, where again the 
Si are Ising spins. The Gaussian distributed Jij have 
mean Jq and variance 1. The SG to FM transition |2^] 
takes place as Jq is increased through the critical value 
Jc ~ 0.96. The mapping from the 2DSG with free bound- 
aries to a general matching problem is given in Ref. [||: 
energy is minimized by selecting a state with a mini- 
mum total weight for frustrated bonds. Timing results 
for the 2DSG are shown in Fig. |3[ One notable difference 
from the RFIM is the apparent convergence to constant 
time per spin in the FM phase. As the algorithm used 
for the 2DSG uses a bond representation, the algorithm 
does not need to distinguish degenerate states related by 
global spin flips. If the frustration is low enough, the FM 
phase is obtained by local operations giving a solution 
with percolating unfrustrated bonds. For low Jo, in the 
SG phase, in contrast, though locally the ground state is 
insensitive to the boundaries, the global ground state is 
sensitive to the disorder and the operations (augmenting 
paths) must therefore be carried out over all scales, so 
that the time per spin may diverge as a power of L (nu- 
merically, approximately as L'^-'^^^^-'^o ^^j. 95 < x < 720 
and Jq — 0.) At the critical point, independent calcu- 
lations [|l^ show that the domain wall fractal dimension 
increases slightly, so that bond updates on large scales 
are more expensive. 

In summary, the time for a polynomial time algorithm 
to find ground states is examined near zero tempera- 
ture critical points for two models, the random field Ising 
model and the 2D spin glass. At the critical points, the 
combinatorial optimization algorithms used, PR fiow and 
general matching, while exponentially faster than, say, 
simulated annealing at finding the ground state, slows 
down. This slowing down of the algorithm is argued to be 
closely related to important physical ideas, namely, the 
uniqueness or two-fold degeneracy of the ground state 
in various phases and the divergence of the correlation 
length as the transition is approached. The time for 
the algorithm can be understood in detail for the RFIM 
with small random fields. It will certainly be of inter- 
est to consider such slowing down and scaling at other 
fixed points for other polynomial time algorithms, to ex- 
pand our understanding of the physics and algorithms for 
these systems. These considerations are clearly relevant 
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for developing efficient parallel algorithms. I would like 
to thank J. Machta for useful discussions. This work was 
supported in part by the NSF (DMR-9702242.) 
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FIG. 1. The CPU time needed to find the ground state 
in the 3D RFIM, for a 766 MHz PHI. The inset plots 
{RL^'^)^^'^S, the scaled number of relabel operations per site, 
with S — {A — Ac)/ Ac, vs. the finite-size scaling variable 
SL^^'^ . The values Ac = 2.27 and v = 1.37 are not fit param- 
eters, but are derived independently. Statistical error bars are 
about 1/5 of the symbol size in both plots. 




FIG. 2. Images of the heights Ui at intermediate stages 
of the PR algorithm, in the 2D RFIM for 100^ samples, for 
A = 10~^ (top) and A = 1.0 (bottom). In the latter case, 
coarsening is cutoff by the finite correlation length; the visible 
spin-up domain shows some structure. The shading is darkest 
at the maximum Ui and is white where Ui — 0. The number 
of global updates, g, is given by the labels. 
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FIG. 3. The CPU time needed to find the ground state in 
the 2D spin glass, as a function of the ferromagnetic strength 
Jo and system size L. 
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